library(AalenJohansen)
library(dplyr)
library(ggplot2)
library(ggpubr)
library(latex2exp)
#### Appendix A.1 ####
M <- t(matrix(c(-3.5,2,1.5,
3,-4,1,
0,0,0),ncol=3,nrow=3))
lambda <- function(t){
1/(1+0.5*t)
}
d_Lambda <- function(t,M,lambda){
lambda(t)*M
}
jump_rate <- function(i,t,u){
#the variable u is not used, it will be used later for time sojourned in the current state
d_L_t <- d_Lambda(t,M,lambda)
J <- dim(d_L_t)[1]
vec <- (1:J==i)*1
-vec%*%d_L_t%*%vec
}
mark_dist <- function(i,s,v){
#the variable v is not used
d_L_t <- d_Lambda(s,M,lambda)
J <- dim(d_L_t)[1]
vec <- (1:J==i)*1
tmp <- (vec %*% d_L_t)* (1-vec)
tmp / sum(tmp)
}
#Simulate paths
simulate_markov_inhomogenous <- function(L) {
R <- runif(L,0,10)
paths <- lapply(1:L,function(n) {
sim_path(1,rates = jump_rate, dist = mark_dist,
tn = R[n], bs= c(R[n]*3,R[n]*4,0))})
}
L <- 10000
set.seed(1)
paths <- simulate_markov_inhomogenous(L)
We can see an example of a path here.
paths[[1]]
## $times
## [1] 0.0000000 0.3246077 0.7910586 2.0599410
##
## $states
## [1] 1 2 1 3
#### Appendix A.2.1 ####
#Convert path data to main_df
paths_to_df <- function(paths){
L <- length(paths)
times <- unlist(lapply(1:L, function(i){
paths[[i]]$times
}))
states <- unlist(lapply(1:L, function(i){
paths[[i]]$states
}))
obs <- unlist(lapply(1:L, function(i){
rep(i,length(paths[[i]]$times))
}))
df <- data.frame(Id = obs, Start_Time = times, Start_State = states)
#End time & end state
df <- df %>%
arrange(Start_Time) %>%
group_by(Id) %>%
mutate(End_Time = data.table::shift(Start_Time,-1),
End_State = data.table::shift(Start_State,-1)) %>%
ungroup() %>%
replace(is.na(.), Inf) %>%
arrange(Id, Start_Time) %>%
mutate(End_State = ifelse(is.infinite(End_State),Start_State,End_State),
Censored = ifelse(Start_State == End_State,
ifelse(is.finite(End_Time),TRUE,FALSE),FALSE)) %>%
group_by(Id) %>%
mutate(Censored = ifelse((!Censored) & is.infinite(End_Time) & (cumsum(Censored) >0),TRUE,Censored)) %>%
ungroup()
df <- df[,c("Id","Start_Time","End_Time","Start_State","End_State","Censored")]
return(df)
}
main_df <- paths_to_df(paths)
We can print the first 10 rows.
main_df[1:10,]
## # A tibble: 10 × 6
## Id Start_Time End_Time Start_State End_State Censored
## <int> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 1 0 0.325 1 2 FALSE
## 2 1 0.325 0.791 2 1 FALSE
## 3 1 0.791 2.06 1 3 FALSE
## 4 1 2.06 Inf 3 3 FALSE
## 5 2 0 0.0591 1 3 FALSE
## 6 2 0.0591 Inf 3 3 FALSE
## 7 3 0 0.128 1 3 FALSE
## 8 3 0.128 Inf 3 3 FALSE
## 9 4 0 0.0784 1 3 FALSE
## 10 4 0.0784 Inf 3 3 FALSE
#### Appendix A.2.2 ####
#Convert main_df to I
df_to_I <- function(df,num_states) {
Init <- df %>%
group_by(State = Start_State) %>%
summarise(Time = 0,
Change = sum(Start_Time == 0))
I <- suppressMessages(df %>%
filter(End_Time < Inf) %>%
mutate(Count_from = TRUE) %>%
bind_rows(filter(.,Censored == FALSE) %>% mutate(Count_from = FALSE)) %>%
arrange(Id,End_Time) %>%
mutate(State = ifelse(Count_from, Start_State,End_State)) %>%
group_by(Time = End_Time, State) %>%
summarise(Change = sum((!Count_from)*1)-sum(Count_from*1)) %>%
ungroup() %>%
bind_rows(Init) %>%
arrange(Time) %>%
group_by(State) %>%
mutate(I_j = cumsum(Change)) %>%
reshape2::dcast(Time ~ State, value.var = "I_j") %>%
zoo::na.locf(na.rm = FALSE) %>%
replace(is.na(.), 0))
if ( sum(!(1:num_states %in% colnames(I))) >0) {
I[,1:num_states[!(1:num_states %in% colnames(I))]] <- 0
}
I <- I[,c("Time",1:num_states)]
I_left_limit <- I
I_left_limit[2:(dim(I)[1]-1),2:dim(I)[2]] <- I[3:dim(I)[1],2:dim(I)[2]]
I_left_limit[dim(I)[1],2:dim(I)[2]] <- I[(dim(I)[1]-1),2:dim(I)[2]]
return(list(I = I, I_left_limit = I_left_limit))
}
I_list <- df_to_I(main_df,3)
The first 10 rows.
I_list$I[1:10,]
## Time 1 2 3
## 1 0.000000e+00 10000 0 0
## 2 8.216051e-05 9999 0 1
## 3 1.200648e-04 9998 1 1
## 4 1.626448e-04 9997 1 2
## 5 3.432192e-04 9996 1 3
## 6 3.554389e-04 9995 2 3
## 7 3.875012e-04 9994 2 4
## 8 3.958256e-04 9993 2 5
## 9 5.169843e-04 9992 2 6
## 10 5.433553e-04 9991 2 7
I_list$I_left_limit[1:10,]
## Time 1 2 3
## 1 0.000000e+00 10000 0 0
## 2 8.216051e-05 9998 1 1
## 3 1.200648e-04 9997 1 2
## 4 1.626448e-04 9996 1 3
## 5 3.432192e-04 9995 2 3
## 6 3.554389e-04 9994 2 4
## 7 3.875012e-04 9993 2 5
## 8 3.958256e-04 9992 2 6
## 9 5.169843e-04 9991 2 7
## 10 5.433553e-04 9990 2 8
#### Appendix A.2.3 ####
#Convert main_df to N
df_to_N <- function(df, num_states) {
N <- df %>%
filter((Censored == FALSE) & (End_Time < Inf)) %>%
arrange(End_Time) %>%
group_by(Start_State, End_State) %>%
mutate(Transitions = cumsum(End_State >= 0),
#j_k = paste0(Start_State,",", End_State)
) %>%
select(Time = End_Time,Start_State, End_State,Transitions) %>%
bind_rows(filter(., Transitions == 1) %>% mutate(Time = 0, Transitions = 0)) %>%
reshape2::dcast(Time ~ Start_State * End_State, value.var = "Transitions") %>%
zoo::na.locf(na.rm = FALSE) %>%
replace(is.na(.), 0)
c <- paste0(unlist(lapply( 1:num_states, function(i) rep(i,num_states))),"_",rep(1:num_states,num_states))
if ( sum(!(c %in% colnames(N))) >0) {
N[,c[!(c %in% colnames(N))]] <- 0
}
N <- N[,c("Time",c)]
for (i in 1:num_states) {
target <- paste0(i,"_",(1:num_states)[!(1:num_states %in% i)])
N[,paste0(i,"_",i)] <- -rowSums(N[,target])
}
delta_N <- N
delta_N[2:dim(N)[1],2:dim(N)[2]] <- delta_N[2:dim(N)[1],2:dim(N)[2]] - delta_N[1:(dim(N)[1]-1),2:dim(N)[2]]
return(list(N = N, delta_N = delta_N))
}
N_list <- df_to_N(main_df, 3)
An example.
N_list$N[1:10,]
## Time 1_1 1_2 1_3 2_1 2_2 2_3 3_1 3_2 3_3
## 1 0.000000e+00 0 0 0 0 0 0 0 0 0
## 2 8.216051e-05 -1 0 1 0 0 0 0 0 0
## 3 1.200648e-04 -2 1 1 0 0 0 0 0 0
## 4 1.626448e-04 -3 1 2 0 0 0 0 0 0
## 5 3.432192e-04 -4 1 3 0 0 0 0 0 0
## 6 3.554389e-04 -5 2 3 0 0 0 0 0 0
## 7 3.875012e-04 -6 2 4 0 0 0 0 0 0
## 8 3.958256e-04 -7 2 5 0 0 0 0 0 0
## 9 5.169843e-04 -8 2 6 0 0 0 0 0 0
## 10 5.433553e-04 -9 2 7 0 0 0 0 0 0
N_list$delta_N[1:10,]
## Time 1_1 1_2 1_3 2_1 2_2 2_3 3_1 3_2 3_3
## 1 0.000000e+00 0 0 0 0 0 0 0 0 0
## 2 8.216051e-05 -1 0 1 0 0 0 0 0 0
## 3 1.200648e-04 -1 1 0 0 0 0 0 0 0
## 4 1.626448e-04 -1 0 1 0 0 0 0 0 0
## 5 3.432192e-04 -1 0 1 0 0 0 0 0 0
## 6 3.554389e-04 -1 1 0 0 0 0 0 0 0
## 7 3.875012e-04 -1 0 1 0 0 0 0 0 0
## 8 3.958256e-04 -1 0 1 0 0 0 0 0 0
## 9 5.169843e-04 -1 0 1 0 0 0 0 0 0
## 10 5.433553e-04 -1 0 1 0 0 0 0 0 0
#### Appendix A.2.4 ####
#Calculate Nelson-Aalen
N_I_to_NA <- function(I_list,N_list, num_states) {
delta_N <- N_list$delta_N
I_left_limit <- I_list$I_left_limit
I_left_limit <- I_left_limit[I_left_limit$Time %in% delta_N$Time,]
I_factor <- delta_N
for (i in 1:num_states) {
target <- paste0(i,"_",1:num_states)
vec <- 1/I_left_limit[,colnames(I_left_limit) == i]
vec <- ifelse(is.infinite(vec),0,vec)
I_factor[,target] <- vec
}
delta_NelsonAalen <- delta_N
delta_NelsonAalen[,2:dim(delta_N)[2]] <- delta_NelsonAalen[,2:dim(delta_N)[2]]*I_factor[,2:dim(delta_N)[2]]
NelsonAalen <- delta_NelsonAalen
NelsonAalen[,2:dim(delta_N)[2]] <- cumsum(delta_NelsonAalen[,2:dim(delta_N)[2]])
return(list(NelsonAalen = NelsonAalen,delta_NelsonAalen=delta_NelsonAalen))
}
NelsonAalen_list <- N_I_to_NA(I_list,N_list,3)
An example.
NelsonAalen_list$NelsonAalen[1:10,]
## Time 1_1 1_2 1_3 2_1 2_2 2_3 3_1 3_2 3_3
## 1 0.000000e+00 0.0000000000 0.00000000 0.0000000000 0 0 0 0 0 0
## 2 8.216051e-05 -0.0001000200 0.00000000 0.0001000200 0 0 0 0 0 0
## 3 1.200648e-04 -0.0002000500 0.00010003 0.0001000200 0 0 0 0 0 0
## 4 1.626448e-04 -0.0003000900 0.00010003 0.0002000600 0 0 0 0 0 0
## 5 3.432192e-04 -0.0004001401 0.00010003 0.0003001100 0 0 0 0 0 0
## 6 3.554389e-04 -0.0005002001 0.00020009 0.0003001100 0 0 0 0 0 0
## 7 3.875012e-04 -0.0006002701 0.00020009 0.0004001801 0 0 0 0 0 0
## 8 3.958256e-04 -0.0007003502 0.00020009 0.0005002602 0 0 0 0 0 0
## 9 5.169843e-04 -0.0008004403 0.00020009 0.0006003502 0 0 0 0 0 0
## 10 5.433553e-04 -0.0009005404 0.00020009 0.0007004503 0 0 0 0 0 0
NelsonAalen_list$delta_NelsonAalen[1:10,]
## Time 1_1 1_2 1_3 2_1 2_2 2_3 3_1 3_2 3_3
## 1 0.000000e+00 0.0000000000 0.00000000 0.0000000000 0 0 0 0 0 0
## 2 8.216051e-05 -0.0001000200 0.00000000 0.0001000200 0 0 0 0 0 0
## 3 1.200648e-04 -0.0001000300 0.00010003 0.0000000000 0 0 0 0 0 0
## 4 1.626448e-04 -0.0001000400 0.00000000 0.0001000400 0 0 0 0 0 0
## 5 3.432192e-04 -0.0001000500 0.00000000 0.0001000500 0 0 0 0 0 0
## 6 3.554389e-04 -0.0001000600 0.00010006 0.0000000000 0 0 0 0 0 0
## 7 3.875012e-04 -0.0001000700 0.00000000 0.0001000700 0 0 0 0 0 0
## 8 3.958256e-04 -0.0001000801 0.00000000 0.0001000801 0 0 0 0 0 0
## 9 5.169843e-04 -0.0001000901 0.00000000 0.0001000901 0 0 0 0 0 0
## 10 5.433553e-04 -0.0001001001 0.00000000 0.0001001001 0 0 0 0 0 0
#### Appendix A.2.5 ####
#Calculate product integral
NA_to_p <- function(I_list,N_list,NelsonAalen_list, num_states) {
#This may be slow
identity <- diag(as.numeric(1), ncol = num_states,nrow=num_states)
Nelson <- NelsonAalen_list$delta_NelsonAalen
Nelson <- lapply(1:dim(Nelson)[1], function(i) matrix(Nelson[i,2:dim(Nelson)[2]],nrow=num_states,ncol = num_states,byrow=TRUE))
Prod_int <- list()
Prod_int[[1]] <- identity
for (i in 2:length(Nelson)) {
Prod_int[[i]] <- Prod_int[[i-1]] %*% (identity + as.numeric(Nelson[[i]]))
}
P <- NelsonAalen_list$delta_NelsonAalen
P[,2:dim(P)[2]] <- matrix(unlist(Prod_int),ncol = num_states**2, byrow = TRUE)
return(P)
}
transitionprobs <- NA_to_p(I_list,N_list,NelsonAalen_list, 3)
An example.
transitionprobs[1:10,]
## Time 1_1 1_2 1_3 2_1 2_2 2_3 3_1 3_2 3_3
## 1 0.000000e+00 1.0000000 0 0 0.00000000 1 0 0.00000000 0 1
## 2 8.216051e-05 0.9999000 0 0 0.00000000 1 0 0.00010002 0 1
## 3 1.200648e-04 0.9998000 0 0 0.00010002 1 0 0.00010002 0 1
## 4 1.626448e-04 0.9996999 0 0 0.00010002 1 0 0.00020004 0 1
## 5 3.432192e-04 0.9995999 0 0 0.00010002 1 0 0.00030006 0 1
## 6 3.554389e-04 0.9994999 0 0 0.00020004 1 0 0.00030006 0 1
## 7 3.875012e-04 0.9993999 0 0 0.00020004 1 0 0.00040008 0 1
## 8 3.958256e-04 0.9992999 0 0 0.00020004 1 0 0.00050010 0 1
## 9 5.169843e-04 0.9991998 0 0 0.00020004 1 0 0.00060012 0 1
## 10 5.433553e-04 0.9990998 0 0 0.00020004 1 0 0.00070014 0 1
#### Appendix A.2.6 ####
#Conditional probabilities
P_conditioned <- function(NelsonAalen_list,s,j,num_states) {
Init <- (1:num_states==j)*1
Nelson <- NelsonAalen_list$delta_NelsonAalen %>%
filter(Time >= s)
Nelson <- lapply(1:dim(Nelson)[1], function(i) matrix(Nelson[i,2:dim(Nelson)[2]],nrow=num_states,ncol = num_states,byrow=TRUE))
Prod_int <- list()
Prod_int[[1]] <- Init
identity <- diag(as.numeric(1), ncol = num_states,nrow=num_states)
for (i in 1:length(Nelson)) {
Prod_int[[i+1]] <- Prod_int[[i]] %*% (identity + as.numeric(Nelson[[i]]))
}
p_con <- NelsonAalen_list$delta_NelsonAalen %>%
filter(Time >= s) %>%
bind_rows(filter(., row_number()==1) %>% mutate(Time = s)) %>%
arrange(Time)
p_con <- p_con[,1:(num_states + 1)]
colnames(p_con) <- c("Time",1:num_states)
p_con[,2:dim(p_con)[2]] <- matrix(unlist(Prod_int),ncol = num_states, byrow = TRUE)
return(p_con)
}
p_con <- P_conditioned(NelsonAalen_list,1,1,3)
An example.
p_con[1:10,]
## Time 1 2 3
## 1 1.000000 1.0000000 0.000000000 0.000000000
## 2 1.000114 1.0000000 0.000000000 0.000000000
## 3 1.000152 0.9994308 0.000000000 0.000569152
## 4 1.000186 0.9988617 0.000569152 0.000569152
## 5 1.000258 0.9982925 0.001138304 0.000569152
## 6 1.000515 0.9977240 0.001138304 0.001137656
## 7 1.000625 0.9977249 0.001137400 0.001137656
## 8 1.000940 0.9971568 0.001705580 0.001137656
## 9 1.001343 0.9971581 0.001704223 0.001137656
## 10 1.001381 0.9971595 0.001702867 0.001137656
#### Appendix A.2.7 ####
#Putting it all together with as-if variant
Estimate <- function(paths,num_states,s= NA, j = NA, as_if = FALSE,debug = TRUE) {
start_time <- Sys.time()
# 1. Start by converting to data frame
if (is.na(s)) {
s <- 0
j <- 1
}
main_df_tmp <- paths_to_df(paths)
if(debug) {
print(paste0("Convert paths to main_df: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
if (as_if) {
tmp <- main_df_tmp %>%
filter(Start_Time <= s) %>%
group_by(Id) %>%
mutate(max_Time = max(Start_Time)) %>%
filter(max_Time == Start_Time) %>%
filter(((is.finite(End_Time) | (Censored == FALSE)))) %>%
filter(Start_State == j)
main_df_tmp <- main_df_tmp %>% filter(Id %in% tmp$Id)
}
# 2. Calculate I and N
I_list_tmp <- df_to_I(main_df_tmp,num_states)
if (debug) {
print(paste0("Convert main_df to I: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
N_list_tmp <- df_to_N(main_df_tmp,num_states)
if (debug) {
print(paste0("Convert main_df to N: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
# 3. Calculate Nelson-Aalen
NelsonAalen_list_tmp <- N_I_to_NA(I_list_tmp,N_list_tmp,num_states)
if (debug) {
print(paste0("Convert I and N to Nelson-Aalen: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
# 4. Calculate Aalen-Johansen
P_tmp <- NA_to_p(I_list_tmp,N_list_tmp,NelsonAalen_list_tmp, 3)
if (debug) {
print(paste0("Convert Nelson-Aalen to transition probs: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
# 5. Calculate probabilties
p_con_tmp <- P_conditioned(NelsonAalen_list_tmp,s,j,num_states)
if (debug) {
print(paste0("Convert transition probs to conditioned probs: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
return(list(I = I_list_tmp$I, I_left_limit = I_list_tmp$I_left_limit,
N = N_list_tmp$N,delta_N = N_list_tmp$delta_N,
NelsonAalen = NelsonAalen_list_tmp$NelsonAalen, delta_NelsonAalen= NelsonAalen_list_tmp$delta_NelsonAalen,
P = P_tmp, p_con = p_con_tmp))
}
total <- Estimate(paths,3)
## [1] "Convert paths to main_df: 0.201 seconds."
## [1] "Convert main_df to I: 0.227 seconds."
## [1] "Convert main_df to N: 0.102 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.011 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.997 seconds."
## [1] "Convert transition probs to conditioned probs: 1.08 seconds."
total <- Estimate(paths,3,s=1,j=1,as_if = TRUE)
## [1] "Convert paths to main_df: 0.198 seconds."
## [1] "Convert main_df to I: 0.135 seconds."
## [1] "Convert main_df to N: 0.023 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.003 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.373 seconds."
## [1] "Convert transition probs to conditioned probs: 0.155 seconds."
#### Appendix A.3 ####
num_states <- 3
s <- 1
j <- 1
debug <- TRUE
pi <- 1
T <- 3
estimate_to_cashflow <- function(estimate,pi,T,n_steps=10000) {
p_con <- estimate$p_con
s <- min(p_con$Time)
t_0 <- p_con[-dim(p_con)[1],1]
t_1 <- p_con[-1,1]
t <- p_con$Time
p_con <- p_con[,2:dim(p_con)[2]]
Nelson <- estimate$delta_NelsonAalen %>%
filter(Time >= s)
A1_increments <- p_con[-1,1]*ifelse(t_1>T,t_1-pmax(T,t_0),0) #A_1(t)-A_1(s)=p(Z_t=1)*(t-s)
A1 <- c(0,cumsum(A1_increments))
A2_increments <- - pi * p_con[-1,1] * ifelse(t_0<T,pmin(T,t_1)-t_0,0) #A_2(t)-A_2(s)=-pi*p(Z_t=1)*(t-s)
A2 <- c(0,cumsum(A2_increments))
A3_increments <- p_con[-1,2] * (t_1-t_0)
A3 <- c(0,cumsum(A3_increments))
A4_increments <- p_con[-1,1] * Nelson[,"1_3"] + p_con[-1,2] * Nelson[,"2_3"]
A4 <- c(0,cumsum(A4_increments))
A <- data.frame(
Time = t,
A1 = A1, A2 = A2, A3 = A3, A4 = A4
) %>% mutate(A = A1 + A2 + A3 + A4)
return(A)
}
estimate <- Estimate(paths,3)
## [1] "Convert paths to main_df: 0.19 seconds."
## [1] "Convert main_df to I: 0.135 seconds."
## [1] "Convert main_df to N: 0.045 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.01 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 1.028 seconds."
## [1] "Convert transition probs to conditioned probs: 1.069 seconds."
cashflow <- estimate_to_cashflow(estimate,1,3)
#Runge-Kutta for true probabilities
runge_kutta <- function(f,a,b,y0,n) {
y <- list()
y[[1]] <- y0
t <- list()
t[[1]] <- a
h <- (b-a)/n
t0 <- a
for (i in 1:n) {
k1 <- f(t0,y0)
k2 <- f(t0 + h/2, y0 + (h/2)*k1)
k3 <- f(t0 + h/2, y0 + (h/2)*k2)
k4 <- f(t0 + h, y0 + h*k3)
y1 <- y0+ (h/6)*(k1+2*k2+2*k3+k4)
y[[i+1]] <- y1
t[[i+1]] <- t0+h
t0 <- t0+h
y0 <- y1
}
return(list(t=t,y=y))
}
n_steps <- 10000
plot_function1 <- function(paths,pi,T,num_states,s= 0, j = 1, debug = TRUE,n_steps=10000) {
#Markov estimate
estimate1 <- Estimate(paths,num_states,s= s, j = j, as_if = FALSE, debug = debug)
start_time <- Sys.time()
cashflow1 <- estimate_to_cashflow(estimate1,pi,T)
if (debug) {
print(paste0("Calculate cashflows: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
#As-If-Markov estimate
estimate2 <- Estimate(paths,num_states,s= s, j = j, as_if = TRUE, debug = debug)
start_time <- Sys.time()
cashflow2 <- estimate_to_cashflow(estimate2,pi,T)
if (debug) {
print(paste0("Calculate cashflows: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
#Plots of probabilities
plotdf <- estimate1$p_con
colnames(plotdf)[2:4] <- paste0("j=",1:3,", Markov")
plotdf <- plotdf%>% reshape2::melt(., id = "Time")
plotdf2 <- estimate2$p_con
colnames(plotdf2)[2:4] <- paste0("j=",1:3,", As-If")
plotdf2 <- plotdf2%>% reshape2::melt(., id = "Time")
#True values (approximated with 4th order runge kutta)
derivative <- function(t,y) {
y %*% d_Lambda(t,M,lambda)
}
y <- runge_kutta(derivative, s,10,(1:num_states == j)*1,n_steps)
plotdf3 <- data.frame(Time = unlist(y$t),
matrix(unlist(y$y),ncol=3,byrow=TRUE))
colnames(plotdf3)[2:4] <- paste0("j=",1:3,", True")
plotdf3 <- plotdf3 %>% reshape2::melt(., id = "Time")
asif_n <- sum(estimate2$I[1,])-estimate2$I[1,1]
markov_n <- sum(estimate1$I[1,])-estimate1$I[1,1]
p1 <- ggplot() +
geom_step(data = plotdf ,mapping = aes(x=Time, y = value,col = variable)) +
geom_step(data = plotdf2,mapping = aes(x=Time, y = value,col = variable),linetype = "dashed") +
geom_line(data = plotdf3,mapping = aes(x=Time, y = value,col = variable),linetype = "dotted",linewidth=1) +
xlim(0,10) +
theme_bw() +
labs(title = TeX(paste0("Occupation probabilities under assumption $\\psi(Z_",s,")=",j,"$ (G)")),
y = TeX(paste0("$P(Z_t=j|G)$")),
x = "t",
subtitle = paste0("As-if estimate based on ",asif_n," observations,\nMarkov based on ",markov_n," observations")) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic")) +
scale_color_manual(values=c("#FA8072", "#FF0000","#8B0000",
"#7B68EE", "#1E90FF","#00008B",
"#3CB371", "#32CD32","#006400"),
name = "Probability, Estimate",
labels = c(TeX("$P(Z_t=a|\\cdot)$, As-If"),
TeX("$P(Z_t=a|\\cdot)$, Markov"),
TeX("$P(Z_t=a|\\cdot)$, True"),
TeX("$P(Z_t=b|\\cdot)$, As-If"),
TeX("$P(Z_t=b|\\cdot)$, Markov"),
TeX("$P(Z_t=b|\\cdot)$, True"),
TeX("$P(Z_t=c|\\cdot)$, As-If"),
TeX("$P(Z_t=c|\\cdot)$, Markov"),
TeX("$P(Z_t=c|\\cdot)$, True")))
#Plots of intensities
plotdf1 <- estimate1$NelsonAalen %>%
select(Time,
`Lambda(1,2), Markov`=`1_2`,
`Lambda(1,3), Markov`=`1_3`,
`Lambda(2,1), Markov`=`2_1`,
`Lambda(2,3), Markov`=`2_3`) %>%
filter(Time >= s)
plotdf1[,2:dim(plotdf1)[2]] <- plotdf1[,2:dim(plotdf1)[2]] - as.data.frame(matrix(as.numeric(rep(plotdf1[1,2:dim(plotdf1)[2]],dim(plotdf1)[1])),ncol = dim(plotdf1)[2]-1,byrow = TRUE))
plotdf1 <- plotdf1 %>%
reshape2::melt(., id = "Time")
plotdf2 <- estimate2$NelsonAalen %>%
select(Time,
`Lambda(1,2), As-if`=`1_2`,
`Lambda(1,3), As-if`=`1_3`,
`Lambda(2,1), As-if`=`2_1`,
`Lambda(2,3), As-if`=`2_3`) %>%
filter(Time >= s)
plotdf2[,2:dim(plotdf2)[2]] <- plotdf2[,2:dim(plotdf2)[2]] - as.data.frame(matrix(as.numeric(rep(plotdf2[1,2:dim(plotdf2)[2]],dim(plotdf2)[1])),ncol = dim(plotdf2)[2]-1,byrow = TRUE))
plotdf2 <- plotdf2 %>%
reshape2::melt(., id = "Time")
times <- s+0:n_steps*(10-s)/n_steps
plotdf3 <- data.frame(Time = times,
matrix(unlist(lapply(times, function(t) as.numeric(t(2*M*(log(1+0.5*t)-log(1+0.5*s)))))),ncol = num_states**2,byrow = TRUE))
colnames(plotdf3) <- c("Time",unlist(lapply(1:num_states, function(i) paste0(i,"_",1:num_states))))
plotdf3 <- plotdf3 %>%
select(Time,
`Lambda(1,2)`=`1_2`,
`Lambda(1,3)`=`1_3`,
`Lambda(2,1)`=`2_1`,
`Lambda(2,3)`=`2_3`) %>%
reshape2::melt(., id = "Time")
p2 <- ggplot() +
geom_step(data = plotdf1 ,mapping = aes(x=Time, y = value,col = variable)) +
geom_step(data = plotdf2,mapping = aes(x=Time, y = value,col = variable),linetype = "dashed") +
geom_line(data = plotdf3,mapping = aes(x=Time, y = value,col = variable),linetype = "dotted",linewidth=1) +
xlim(0,10) +
theme_bw() +
labs(title = TeX(paste0("Cumulative transition rates under assumption $\\psi(Z_",s,")=",j,"$ (G)")),
y = TeX(paste0("$\\Lambda_{jk}(t|G)-\\Lambda_{jk}(s|G)$")),
x = "t",
subtitle = paste0("As-if estimate based on ",asif_n," observations,\nMarkov based on ",markov_n," observations")) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic")) +
scale_color_manual(values=c("#FA8072", "#FF0000","#8B0000","#7B68EE",
"#1E90FF","#00008B","#3CB371", "#32CD32",
"#006400","#D2691E","#F4A460","#DEB887"),
name = "j to k, Estimate",
labels = c("j=a, k=b, True",
"j=a, k=b, As-If",
"j=a, k=b, Markov",
"j=a, k=c, True",
"j=a, k=c, As-If",
"j=a, k=c, Markov",
"j=b, k=a, True",
"j=b, k=a, As-If",
"j=b, k=a, Markov",
"j=b, k=c, True",
"j=b, k=c, As-If",
"j=b, k=c, Markov"))
#Plots of cashflows
colnames(cashflow1)[2:dim(cashflow1)[2]] <- paste0(colnames(cashflow1)[2:dim(cashflow1)[2]],", Markov")
plotdf1 <- cashflow1 %>% reshape2::melt(id = "Time")
colnames(cashflow2)[2:dim(cashflow2)[2]] <- paste0(colnames(cashflow2)[2:dim(cashflow2)[2]],", As-If")
plotdf2 <- cashflow2 %>% reshape2::melt(id = "Time")
#Calculate true values
prob_to_1 <- approxfun(y$t,y=unlist(y$y)[1:length(y$t)*length(y$y[[1]])-2])
prob_to_2 <- approxfun(y$t,y=unlist(y$y)[1:length(y$t)*length(y$y[[1]])-1])
prob_to_3 <- approxfun(y$t,y=unlist(y$y)[1:length(y$t)*length(y$y[[1]])-0])
A1_derivative <- function(t,y){
(t>T)*prob_to_1(t)
}
A1 <- runge_kutta(A1_derivative,s,10,0,n_steps)
A2_derivative <- function(t,y){
-pi*(t<=T)*prob_to_1(t)
}
A2 <- runge_kutta(A2_derivative,s,10,0,n_steps)
A3_derivative <- function(t,y){
prob_to_2(t)
}
A3 <- runge_kutta(A3_derivative,s,10,0,n_steps)
A4_derivative <- function(t,y){
#prob_to_1(t-)=prob_to_1(t) same for prob_to_2
prob_to_1(t)*d_Lambda(t,M,lambda)[1,3]+prob_to_2(t)*d_Lambda(t,M,lambda)[2,3]
}
A4 <- runge_kutta(A4_derivative,s,10,0,n_steps)
cashflow3 <- data.frame(
Time = unlist(A1$t),
A1 = unlist(A1$y),
A2 = unlist(A2$y),
A3 = unlist(A3$y),
A4 = unlist(A4$y)
) %>%
mutate(A = A1 + A2 + A3 + A4)
colnames(cashflow3) <- c("Time","A1_true","A2_true","A3_true","A4_true","A_true")
plotdf3 <- cashflow3 %>% reshape2::melt(id = "Time")
if (debug) {
print(paste0("Calculate theoretical cashflows and more: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
p3 <- ggplot() +
geom_step(data = plotdf1 ,mapping = aes(x=Time, y = value,col = variable)) +
geom_step(data = plotdf2,mapping = aes(x=Time, y = value,col = variable), linetype = "dashed") +
geom_line(data = plotdf3,mapping = aes(x=Time, y = value,col = variable), linetype = "dotted", linewidth=1) +
geom_vline(xintercept = T, col = "black",linetype = "dashed") +
xlim(0,10) +
theme_bw() +
labs(title = TeX(paste0("Accumulated cash-flow under assumption $\\psi(Z_",s,")=",j,"$ (G)")),
y = TeX(paste0("A(t|G)-A(s|G)")),
x = "t",
subtitle = paste0("As-if estimate based on ",asif_n," observations,\nMarkov based on ",markov_n," observations")) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic")) +
scale_color_manual(values=c("#FA8072", "#FF0000","#8B0000","#7B68EE",
"#1E90FF","#00008B","#3CB371", "#32CD32",
"#006400","#D2691E","#F4A460","#DEB887",
"#FF1493","#C71585","#FFB6C1"),
name = "Cash-flow, Estimate",
labels = c(TeX("Total ($A$), Markov"),
TeX("Total ($A$), As-If"),
TeX("Total ($A$), True"),
TeX("Pension ($A_1$), Markov"),
TeX("Pension ($A_1$), As-If"),
TeX("Pension ($A_1$), True"),
TeX("Premium ($A_2$), Markov"),
TeX("Premium ($A_2$), As-If"),
TeX("Premium ($A_2$), True"),
TeX("Disabled ($A_3$), Markov"),
TeX("Disabled ($A_3$), As-If"),
TeX("Disabled ($A_3$), True"),
TeX("Death ($A_4$), Markov"),
TeX("Death ($A_4$), As-If"),
TeX("Death ($A_4$), True")))
return(list(p1= p1, p2 = p2, p3 = p3))
}
plot1 <- plot_function1(paths,1,3,3,s= 0, j = 1)
## [1] "Convert paths to main_df: 0.324 seconds."
## [1] "Convert main_df to I: 0.128 seconds."
## [1] "Convert main_df to N: 0.044 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.03 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.98 seconds."
## [1] "Convert transition probs to conditioned probs: 1.135 seconds."
## [1] "Calculate cashflows: 0.003 seconds."
## [1] "Convert paths to main_df: 0.193 seconds."
## [1] "Convert main_df to I: 0.243 seconds."
## [1] "Convert main_df to N: 0.046 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.009 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 1.061 seconds."
## [1] "Convert transition probs to conditioned probs: 1.028 seconds."
## [1] "Calculate cashflows: 0.003 seconds."
## [1] "Calculate theoretical cashflows and more: 0.783 seconds."
plot2 <- plot_function1(paths,1,3,3,s= 1, j = 1)
## [1] "Convert paths to main_df: 0.202 seconds."
## [1] "Convert main_df to I: 0.127 seconds."
## [1] "Convert main_df to N: 0.044 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.024 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 1.054 seconds."
## [1] "Convert transition probs to conditioned probs: 0.296 seconds."
## [1] "Calculate cashflows: 0.002 seconds."
## [1] "Convert paths to main_df: 0.236 seconds."
## [1] "Convert main_df to I: 0.152 seconds."
## [1] "Convert main_df to N: 0.023 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.003 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.326 seconds."
## [1] "Convert transition probs to conditioned probs: 0.196 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
## [1] "Calculate theoretical cashflows and more: 0.684 seconds."
plot3 <- plot_function1(paths,1,3,3,s= 3, j = 1)
## [1] "Convert paths to main_df: 0.241 seconds."
## [1] "Convert main_df to I: 0.252 seconds."
## [1] "Convert main_df to N: 0.042 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.007 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 1.004 seconds."
## [1] "Convert transition probs to conditioned probs: 0.057 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
## [1] "Convert paths to main_df: 0.309 seconds."
## [1] "Convert main_df to I: 0.111 seconds."
## [1] "Convert main_df to N: 0.017 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.001 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.087 seconds."
## [1] "Convert transition probs to conditioned probs: 0.026 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
## [1] "Calculate theoretical cashflows and more: 0.637 seconds."
plot4 <- plot_function1(paths,1,3,3,s= 6, j = 1)
## [1] "Convert paths to main_df: 0.203 seconds."
## [1] "Convert main_df to I: 0.136 seconds."
## [1] "Convert main_df to N: 0.042 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.007 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 1.057 seconds."
## [1] "Convert transition probs to conditioned probs: 0.007 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
## [1] "Convert paths to main_df: 0.247 seconds."
## [1] "Convert main_df to I: 0.12 seconds."
## [1] "Convert main_df to N: 0.017 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.001 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.019 seconds."
## [1] "Convert transition probs to conditioned probs: 0.004 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
## [1] "Calculate theoretical cashflows and more: 0.755 seconds."
p1 <- ggarrange(plotlist = list(plot1$p1,plot2$p1,plot3$p1,plot4$p1),ncol = 2,nrow=2)
## Warning: Removed 3 rows containing missing values (`geom_line()`).
## Removed 3 rows containing missing values (`geom_line()`).
p2 <- ggarrange(plotlist = list(plot1$p2,plot2$p2,plot3$p2,plot4$p2),ncol = 2,nrow=2)
p3 <- ggarrange(plotlist = list(plot1$p3,plot2$p3,plot3$p3,plot4$p3),ncol = 2,nrow=2)
## Warning: Removed 5 rows containing missing values (`geom_line()`).
## Warning: Removed 5 rows containing missing values (`geom_line()`).
scaler <- 1.75
ggsave("plotA1.png",p1,units = "px", width = 1920*scaler,height = 1080*scaler,scale = 1.5)
ggsave("plotA2.png",p2,units = "px", width = 1920*scaler,height = 1080*scaler,scale = 1.5)
ggsave("plotA3.png",p3,units = "px", width = 1920*scaler,height = 1080*scaler,scale = 1.5)
We can show the plots below.
knitr::include_graphics("plotA1.png")
knitr::include_graphics("plotA2.png")
knitr::include_graphics("plotA3.png")
#### Appendix A.4 ####
L <- ceiling(2*65**(2*0:100/100+1))
S <- 6
J <- 1
supremums <- function(L,S,J, debug = TRUE,n_steps = 10000) {
results <- data.frame(L = NA,s=NA,j=NA,as_if = TRUE,`Lambda(1,2)`=NA,`Lambda(1,3)`=NA,`Lambda(2,1)`=NA,`Lambda(2,3)`=NA)
L <- as.integer(L)
counter <- 1
set.seed(1)
if (debug) {
print(paste0(Sys.time()," Starting simulating ",max(L)," sample paths."))
}
paths <- simulate_markov_inhomogenous(max(L))
if (debug) {
print(paste0(Sys.time()," Done simulating ",max(L)," sample paths."))
}
for (l in L) {
tmp_paths <- paths[1:l]
for (s in S) {
for (j in J) {
#Generate paths
estimate1 <- tryCatch(Estimate(tmp_paths,3, s= s, j=j,as_if = FALSE, debug = FALSE),
error = function(e) NA)
estimate2 <- tryCatch(Estimate(tmp_paths,3,s= s, j=j,as_if = TRUE, debug = FALSE),
error = function(e) NA)
if ((length(estimate1)>1) & (length(estimate2)>1)) {
markov <- estimate1$NelsonAalen %>%
select(Time,
`Lambda(1,2)`=`1_2`,
`Lambda(1,3)`=`1_3`,
`Lambda(2,1)`=`2_1`,
`Lambda(2,3)`=`2_3`) %>%
filter(Time >= s)
markov[,2:dim(markov)[2]] <- markov[,2:dim(markov)[2]] -
as.data.frame(matrix(as.numeric(rep(markov[1,2:dim(markov)[2]],dim(markov)[1])),ncol = dim(markov)[2]-1,byrow = TRUE))
as_if_markov <- estimate2$NelsonAalen %>%
select(Time,
`Lambda(1,2)`=`1_2`,
`Lambda(1,3)`=`1_3`,
`Lambda(2,1)`=`2_1`,
`Lambda(2,3)`=`2_3`) %>%
filter(Time >= s)
as_if_markov[,2:dim(as_if_markov)[2]] <- as_if_markov[,2:dim(as_if_markov)[2]] -
as.data.frame(matrix(as.numeric(rep(as_if_markov[1,2:dim(as_if_markov)[2]],dim(as_if_markov)[1])),ncol = dim(as_if_markov)[2]-1,byrow = TRUE))
times <- sort(unique(c(s+0:n_steps*(10-s)/n_steps,as_if_markov$Time,markov$Time)),decreasing = FALSE)
true_values <- data.frame(Time = times,
matrix(unlist(lapply(times, function(t) as.numeric(t(2*M*(log(1+0.5*t)-log(1+0.5*s)))))),ncol = num_states**2,byrow = TRUE))
colnames(true_values) <- c("Time",unlist(lapply(1:num_states, function(i) paste0(i,"_",1:num_states))))
true_values <- true_values %>%
select(Time,
`Lambda(1,2),true`=`1_2`,
`Lambda(1,3),true`=`1_3`,
`Lambda(2,1),true`=`2_1`,
`Lambda(2,3),true`=`2_3`)
markov <- merge(markov,true_values,all.x = TRUE)
as_if_markov <- merge(as_if_markov,true_values,all.x = TRUE)
results[counter,] <- c(
l,s,j,FALSE,
max(abs(markov$`Lambda(1,2)`-markov$`Lambda(1,2),true`)),
max(abs(markov$`Lambda(1,3)`-markov$`Lambda(1,3),true`)),
max(abs(markov$`Lambda(2,1)`-markov$`Lambda(2,1),true`)),
max(abs(markov$`Lambda(2,3)`-markov$`Lambda(2,3),true`))
)
results[counter+1,] <- c(
l,s,j,TRUE,
max(abs(as_if_markov$`Lambda(1,2)`-as_if_markov$`Lambda(1,2),true`)),
max(abs(as_if_markov$`Lambda(1,3)`-as_if_markov$`Lambda(1,3),true`)),
max(abs(as_if_markov$`Lambda(2,1)`-as_if_markov$`Lambda(2,1),true`)),
max(abs(as_if_markov$`Lambda(2,3)`-as_if_markov$`Lambda(2,3),true`))
)
counter <- counter + 2
if (debug == TRUE) {
print(paste0(Sys.time(),": Done with L=",l,", s=",s," and j=",j,"."))
}
}
}
}
}
return(results)
}
plot_function2 <- function(results, s, j) {
plotdf <- results %>% filter((s==s) & (j==j)) %>% filter(as_if == TRUE) %>%
select(-s,-j,-as_if) %>%
melt(id = "L")
plotdf2 <- results %>% filter((s==s) & (j==j)) %>% filter(as_if == FALSE) %>%
select(-s,-j,-as_if) %>%
melt(id = "L")
ggplot() + geom_point(data = plotdf %>% filter(variable == "Lambda.1.3."), aes(x=L,y = value, col = variable))+
geom_line(data = plotdf2 %>% filter(variable == "Lambda.1.3."), aes(x=L,y = value, col = variable)) +
scale_x_continuous(trans = "log2") + scale_y_continuous(trans = "log2")
}
#Not run
#results1 <- supremums(L = 1:10*1000,S=0:6,J=1) #run-time ≈
#results2 <- supremums(L = ceil(2*65**(2*0:25/25+1)),S=6,J=1,n_steps = 1000000) #run-time ≈ 11.5 min
#results3 <- supremums(L = ceil(2*88**(2*0:25/25+1)),S=6,J=2) #run-time ≈
#plot1 <- plot_function(results,6,1)
#plot1
#Calcuæate P(Z_s=j)P(R>s)
s <- 6
c(1,0,0)%*%expm::expm(M*2*log(1+0.5*s))
## [,1] [,2] [,3]
## [1,] 0.01549959 0.01142945 0.973071
##################################
###### Estimate Semi-Markov ######
##################################
#### Appendix B.1 ####
lambda_a_b <- function(t,u) {
0.09+0.0012*t+(t>u)*0.015*t
}
lambda_a_c <- function(t,u) {
0.01+0.002*t+(t>u)*0.001*t
}
lambda_b_a <- function(t,u) {
0.05/max(1,u**2)
}
lambda_b_c <- function(t,u) {
0.03+0.002*t+0.025*4/max(4,u)
}
jump_rate <- function(i, t, u){
if(i == 1){
lambda_a_b(t,u)+lambda_a_c(t,u)
} else if(i == 2){
lambda_b_a(t,u)+lambda_b_c(t,u)
} else{
0
}
}
mark_dist <- function(i, s, v){
if(i == 1){
tmp1 <- lambda_a_b(s,v)
tmp2 <- lambda_a_c(s,v)
c(0, tmp1/(tmp1+tmp2), tmp2/(tmp1+tmp2))
} else if(i == 2){
tmp1 <- lambda_b_a(s,v)
tmp2 <- lambda_b_c(s,v)
c(tmp1/(tmp1+tmp2),0, tmp2/(tmp1+tmp2))
} else{
0
}
}
#Simulate paths
simulate_semimarkov <- function(L) {
R <- runif(L,10,40)
paths <- lapply(1:L,function(n) {
sim_path(1,rates = jump_rate, dist = mark_dist,tn = R[n],
bs = c(lambda_a_b(R[n],0)+lambda_a_c(R[n],0),
lambda_b_a(R[n],0)+lambda_b_c(R[n],0),
0))})
}
L <- 10000
set.seed(1)
paths <- simulate_semimarkov(L)
We can see an example of a path here.
paths[[1]]
## $times
## [1] 0.000000 5.811178 17.965260
##
## $states
## [1] 1 2 2
#### Appendix B.2 ####
num_states <- 3
s <- 20
j <- 1
debug <- TRUE
pi <- 1
T <- 15
plot_function_B1 <- function(paths,pi,T,num_states,s= 0, j = 1, debug = TRUE,n_steps=10000) {
#Markov estimate
estimate1 <- Estimate(paths,num_states,s= s, j = j, as_if = FALSE, debug = debug)
start_time <- Sys.time()
cashflow1 <- estimate_to_cashflow(estimate1,pi,T)
if (debug) {
print(paste0("Calculate cashflows: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
#As-If-Markov estimate
estimate2 <- Estimate(paths,num_states,s= s, j = j, as_if = TRUE, debug = debug)
start_time <- Sys.time()
cashflow2 <- estimate_to_cashflow(estimate2,pi,T)
if (debug) {
print(paste0("Calculate cashflows: ",round(as.numeric(Sys.time()-start_time,units = "secs"),digits=3)," seconds."))
start_time <- Sys.time()
}
#Plots of probabilities
plotdf <- estimate1$p_con
colnames(plotdf)[2:4] <- paste0("j=",1:3,", Markov")
plotdf <- plotdf%>% reshape2::melt(., id = "Time")
plotdf2 <- estimate2$p_con
colnames(plotdf2)[2:4] <- paste0("j=",1:3,", As-If")
plotdf2 <- plotdf2%>% reshape2::melt(., id = "Time")
asif_n <- sum(estimate2$I[1,])-estimate2$I[1,1]
markov_n <- sum(estimate1$I[1,])-estimate1$I[1,1]
p1 <- ggplot() +
geom_step(data = plotdf ,mapping = aes(x=Time, y = value,col = variable)) +
geom_step(data = plotdf2,mapping = aes(x=Time, y = value,col = variable),linetype = "dashed") +
#geom_line(data = plotdf3,mapping = aes(x=Time, y = value,col = variable),linetype = "dotted",linewidth=1) +
xlim(0,40) +
theme_bw() +
labs(title = TeX(paste0("Occupation probabilities under assumption $\\psi(Z_{",s,"})=",j,"$ (G)")),
y = TeX(paste0("$P(Z_t=j|G)$")),
x = "t",
subtitle = paste0("As-if estimate based on ",asif_n," observations,\nMarkov based on ",markov_n," observations")) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic")) +
scale_color_manual(values=c("#FA8072", "#FF0000",#"#8B0000",
"#7B68EE", "#1E90FF",#"#00008B",
"#3CB371", "#32CD32"#,"#006400"
),
name = "Probability, Estimate",
labels = c(TeX("$P(Z_t=a|\\cdot)$, As-If"),
TeX("$P(Z_t=a|\\cdot)$, Markov"),
#TeX("$P(Z_t=a|\\cdot)$, True"),
TeX("$P(Z_t=b|\\cdot)$, As-If"),
TeX("$P(Z_t=b|\\cdot)$, Markov"),
#TeX("$P(Z_t=b|\\cdot)$, True"),
TeX("$P(Z_t=c|\\cdot)$, As-If"),
TeX("$P(Z_t=c|\\cdot)$, Markov")
#TeX("$P(Z_t=c|\\cdot)$, True")
))
#Plots of intensities
plotdf1 <- estimate1$NelsonAalen %>%
select(Time,
`Lambda(1,2), Markov`=`1_2`,
`Lambda(1,3), Markov`=`1_3`,
`Lambda(2,1), Markov`=`2_1`,
`Lambda(2,3), Markov`=`2_3`) %>%
filter(Time >= s)
plotdf1[,2:dim(plotdf1)[2]] <- plotdf1[,2:dim(plotdf1)[2]] - as.data.frame(matrix(as.numeric(rep(plotdf1[1,2:dim(plotdf1)[2]],dim(plotdf1)[1])),ncol = dim(plotdf1)[2]-1,byrow = TRUE))
plotdf1 <- plotdf1 %>%
reshape2::melt(., id = "Time")
plotdf2 <- estimate2$NelsonAalen %>%
select(Time,
`Lambda(1,2), As-if`=`1_2`,
`Lambda(1,3), As-if`=`1_3`,
`Lambda(2,1), As-if`=`2_1`,
`Lambda(2,3), As-if`=`2_3`) %>%
filter(Time >= s)
plotdf2[,2:dim(plotdf2)[2]] <- plotdf2[,2:dim(plotdf2)[2]] - as.data.frame(matrix(as.numeric(rep(plotdf2[1,2:dim(plotdf2)[2]],dim(plotdf2)[1])),ncol = dim(plotdf2)[2]-1,byrow = TRUE))
plotdf2 <- plotdf2 %>%
reshape2::melt(., id = "Time")
p2 <- ggplot() +
geom_step(data = plotdf1 ,mapping = aes(x=Time, y = value,col = variable)) +
geom_step(data = plotdf2,mapping = aes(x=Time, y = value,col = variable),linetype = "dashed") +
#geom_line(data = plotdf3,mapping = aes(x=Time, y = value,col = variable),linetype = "dotted",linewidth=1) +
xlim(0,40) +
theme_bw() +
labs(title = TeX(paste0("Cumulative transition rates under assumption $\\psi(Z_{",s,"})=",j,"$ (G)")),
y = TeX(paste0("$\\Lambda_{jk}(t|G)-\\Lambda_{jk}(s|G)$")),
x = "t",
subtitle = paste0("As-if estimate based on ",asif_n," observations,\nMarkov based on ",markov_n," observations")) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic")) +
scale_color_manual(values=c("#FA8072", "#FF0000",#"#8B0000",
"#7B68EE","#1E90FF",#"#00008B",
"#3CB371", "#32CD32",#"#006400",
"#D2691E","#F4A460"#,"#DEB887"
),
name = "j to k, Estimate",
labels = c(#"j=a, k=b, True",
"j=a, k=b, As-If",
"j=a, k=b, Markov",
#"j=a, k=c, True",
"j=a, k=c, As-If",
"j=a, k=c, Markov",
#"j=b, k=a, True",
"j=b, k=a, As-If",
"j=b, k=a, Markov",
#"j=b, k=c, True",
"j=b, k=c, As-If",
"j=b, k=c, Markov"))
#Plots of cashflows
colnames(cashflow1)[2:dim(cashflow1)[2]] <- paste0(colnames(cashflow1)[2:dim(cashflow1)[2]],", Markov")
plotdf1 <- cashflow1 %>% reshape2::melt(id = "Time")
colnames(cashflow2)[2:dim(cashflow2)[2]] <- paste0(colnames(cashflow2)[2:dim(cashflow2)[2]],", As-If")
plotdf2 <- cashflow2 %>% reshape2::melt(id = "Time")
p3 <- ggplot() +
geom_step(data = plotdf1 ,mapping = aes(x=Time, y = value,col = variable)) +
geom_step(data = plotdf2,mapping = aes(x=Time, y = value,col = variable), linetype = "dashed") +
#geom_line(data = plotdf3,mapping = aes(x=Time, y = value,col = variable), linetype = "dotted", linewidth=1) +
geom_vline(xintercept = T, col = "black",linetype = "dashed") +
xlim(0,40) +
theme_bw() +
labs(title = TeX(paste0("Accumulated cash-flow under assumption $\\psi(Z_{",s,"})=",j,"$ (G)")),
y = TeX(paste0("A(t|G)-A(s|G)")),
x = "t",
subtitle = paste0("As-if estimate based on ",asif_n," observations,\nMarkov based on ",markov_n," observations")) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(face = "italic")) +
scale_color_manual(values=c("#FA8072", "#FF0000",#"#8B0000",
"#7B68EE","#1E90FF",#"#00008B",
"#3CB371", "#32CD32",#"#006400",
"#D2691E","#F4A460",#"#DEB887",
"#FF1493","#C71585"#,"#FFB6C1"
),
name = "Cash-flow, Estimate",
labels = c(TeX("Total ($A$), As-If"),
TeX("Total ($A$), Markov"),
TeX("Pension ($A_1$), As-If"),
TeX("Pension ($A_1$), Markov"),
TeX("Premium ($A_2$), As-If"),
TeX("Premium ($A_2$), Markov"),
TeX("Disabled ($A_3$), As-If"),
TeX("Disabled ($A_3$), Markov"),
TeX("Death ($A_4$), As-If"),
TeX("Death ($A_4$), Markov")))
return(list(p1= p1, p2 = p2, p3 = p3))
}
plot1 <- plot_function_B1(paths,pi = 1, T= 15,num_states = 3,s= 0, j = 1)
## [1] "Convert paths to main_df: 0.193 seconds."
## [1] "Convert main_df to I: 0.095 seconds."
## [1] "Convert main_df to N: 0.033 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.006 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.637 seconds."
## [1] "Convert transition probs to conditioned probs: 0.725 seconds."
## [1] "Calculate cashflows: 0.002 seconds."
## [1] "Convert paths to main_df: 0.214 seconds."
## [1] "Convert main_df to I: 0.209 seconds."
## [1] "Convert main_df to N: 0.03 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.005 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.69 seconds."
## [1] "Convert transition probs to conditioned probs: 0.802 seconds."
## [1] "Calculate cashflows: 0.003 seconds."
plot2 <- plot_function_B1(paths,pi = 1, T= 15,num_states = 3,s= 10, j = 1)
## [1] "Convert paths to main_df: 0.22 seconds."
## [1] "Convert main_df to I: 0.105 seconds."
## [1] "Convert main_df to N: 0.039 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.006 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.899 seconds."
## [1] "Convert transition probs to conditioned probs: 0.253 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
## [1] "Convert paths to main_df: 0.187 seconds."
## [1] "Convert main_df to I: 0.124 seconds."
## [1] "Convert main_df to N: 0.019 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.002 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.174 seconds."
## [1] "Convert transition probs to conditioned probs: 0.158 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
plot3 <- plot_function_B1(paths,pi = 1, T= 15,num_states = 3,s= 10, j = 2)
## [1] "Convert paths to main_df: 0.238 seconds."
## [1] "Convert main_df to I: 0.1 seconds."
## [1] "Convert main_df to N: 0.031 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.005 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.688 seconds."
## [1] "Convert transition probs to conditioned probs: 0.286 seconds."
## [1] "Calculate cashflows: 0.002 seconds."
## [1] "Convert paths to main_df: 0.204 seconds."
## [1] "Convert main_df to I: 0.193 seconds."
## [1] "Convert main_df to N: 0.022 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.003 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.314 seconds."
## [1] "Convert transition probs to conditioned probs: 0.115 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
plot4 <- plot_function_B1(paths,pi = 1, T= 15,num_states = 3,s= 20, j = 2)
## [1] "Convert paths to main_df: 0.227 seconds."
## [1] "Convert main_df to I: 0.108 seconds."
## [1] "Convert main_df to N: 0.032 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.005 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.81 seconds."
## [1] "Convert transition probs to conditioned probs: 0.083 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
## [1] "Convert paths to main_df: 0.235 seconds."
## [1] "Convert main_df to I: 0.148 seconds."
## [1] "Convert main_df to N: 0.021 seconds."
## [1] "Convert I and N to Nelson-Aalen: 0.001 seconds."
## [1] "Convert Nelson-Aalen to transition probs: 0.17 seconds."
## [1] "Convert transition probs to conditioned probs: 0.06 seconds."
## [1] "Calculate cashflows: 0.001 seconds."
p1 <- ggarrange(plotlist = list(plot2$p1,plot3$p1,plot2$p2,plot3$p2,plot2$p3,plot3$p3),ncol = 2,nrow=3)
scaler <- 1.75
ggsave("plotB1.png",p1,units = "px", width = 1920*scaler,height = 1920*scaler,scale = 1.5)
We can show the plots below.
knitr::include_graphics("plotB1.png")